Inst/doc/mixed von mises.R

In Atan1988/FlexFit: Fits Severity distributions for Insurance data

library(circular)
library(dplyr)

data("forAllen2")

dmvonmises_unif <- function(x, prop_unif, prop1, prop2, mu1, mu2, kappa1, kappa2) {
 a <- -pi
 b <- pi

 p1 <- exp(prop_unif) / (exp(prop_unif) + exp(prop1) + exp(prop2))
 p2 <- exp(prop1) / (exp(prop_unif) + exp(prop1) + exp(prop2))
 p3 <- exp(prop2) / (exp(prop_unif) + exp(prop1) + exp(prop2))

 prop_mvonmises <- p2 + p3

 p1 * dunif(x, min = a, max = b) +
   prop_mvonmises * circular::dmixedvonmises(x, mu1, mu2, kappa1, kappa2, prop = p2 / prop_mvonmises)
}

dmvonmises_unif2 <- function(x, prop_unif, prop1, prop2, mu1, kappa1) {
  a <- -pi
  b <- pi

  p1 <- exp(prop_unif) / (exp(prop_unif) + exp(prop1) + exp(prop2))
  p2 <- exp(prop1) / (exp(prop_unif) + exp(prop1) + exp(prop2))
  p3 <- exp(prop2) / (exp(prop_unif) + exp(prop1) + exp(prop2))

  prop_mvonmises <- p2 + p3

  p1 * dunif(x, min = a, max = b) +
    prop_mvonmises * circular::dmixedvonmises(x, mu1, mu2 = mu1 + pi, kappa1, kappa1, prop = p2 / prop_mvonmises)
}

start_pt <- c(prop_unif = 0,
              prop1 = .348,
              prop2 = .652,
              mu1 = 0.041,
              mu2 = 0.041 + pi,
              kappa1 = 3.48,
              kappa2 = 3.48
              )

start_pt2 <- c(prop_unif = 0,
              prop1 = .348,
              prop2 = .652,
              mu1 = 0.041,
              kappa1 = 3.48
)

x1 <- forAllen2$X1

system.time(fit1b <- sevfit(x1, Att = -pi, FUN = 'mvonmises_unif',
                            ini = start_pt,
                            parnames = names(start_pt), gn_ll_mtd = NULL,
                            hessian = T
)
)

system.time(fit2b <- sevfit(x1, Att = -pi, FUN = 'mvonmises_unif2',
                            ini = start_pt2,
                            parnames = names(start_pt2), gn_ll_mtd = NULL,
                            hessian = T
)
)

suppressWarnings(
  logLik(fit2b, y = x1)
); suppressWarnings(
  logLik(fit1b, y = x1)
); suppressWarnings(
  coef(fit2b)
); suppressWarnings(
  coef(fit1b)
)

c((fit1b$parameters[1:3] %>% unlist())/ sum(fit1b$parameter[1:3] %>% unlist()),
fit1b$parameters[4:7] %>% unlist())


c((fit2b$parameters[1:3] %>% unlist())/ sum(fit2b$parameter[1:3] %>% unlist()),
  fit2b$parameters[4:7] %>% unlist())


dmvonmises_unif2(forAllen$X1, prop_unif = 0, prop1 = 0.4919,
                 prop2 = 0.5081, mu1 = -0.0729 , kappa1 = 3.0370) %>% log() %>% sum()

suppressWarnings(
  dmvonmises_unif2(forAllen$X1, prop_unif = 0.3512, prop1 = 0.2904,
                 prop2 = 0.3587, mu1 = -0.0733, kappa1 = 6.0973) %>% log() %>% sum()
)
dmvonmises_unif2(forAllen$X1, prop_unif = -22.03645, prop1 = -46.25202,
                 prop2 = -47.7848, mu1 = -0.08692372, kappa1 = 6.781965) %>% log() %>% sum()


suppressWarnings(
  do.call(dmvonmises_unif2, append(list(x = forAllen$X1),  fit2b$parameters)) %>% log() %>% sum()
)